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The relativistic correction to the QCD static inter-quark potential at 0(l/m) is investigated 
nonperturbatively for the first time by using lattice Monte Carlo QCD simulations. The correc- 
tion is found to be comparable with the Coulombic term of the static potential when applied to 
charmonium, and amounts to one-fourth of the Coulombic term for bottomonium. 



PACS numbers: 11.15.Ha, 12.38.Gc, 12.39.Pn 



I. INTRODUCTION 

Heavy quarkonia, i.e. bound states of a heavy quark 
and antiquark 0, H H H, offer a unique opportunity 
to gain an understanding of nonperturbative QCD. A 
possible way of studying such systems systematically in 
QCD is to employ nonrelativistic QCD (NRQCD) [IQ, 
which is obtained by integrating out the scale above the 
heavy quark mass m ^> Aqcd- Further, by integrat- 
ing out the scale mv, where v is quark velocity, one 
arrives at a framework called potential NRQCD (pN- 
RQCD) 0,0 EJ El i where the static potential emerges 
as the leading-order contribution, followed by relativis- 
tic corrections in powers of 1/m. The potential at 
0(l/m 2 ) contains the leading order spin-dependent cor- 
rections Qjl IeH El and the velocity-dependent poten- 
tials fl4lll5l|. Perturbation theory may be applied to the 
determination of these potentials to some extent. How- 
ever, since the binding energy is typically of the scale 
mi) 2 , which can be of the same order as Aqcd due to 
the nonrelativistic nature of the system, v <C 1, as well 
as the fact that perturbation theory cannot incorporate 
quark confinement, it is essential to determine the poten- 
tial nonperturbatively. The various properties of heavy 
quarkonium can be extracted by solving the Schrodingcr 
equation with these potentials. 

Monte Carlo simulations of lattice QCD offer a pow- 
erful tool for the nonperturbative determination of the 
potentials, and it is the aim of this Letter to present the 
simulation result of the heavy quark potential at 0(l/m), 
which has not been investigated so far on the lattice. Let 
us denote the spatial position of the quark and antiquark 
as ri and with the relative distance r — \r\ — f%\ and 
the masses mi and m-i, respectively. The potential is 

V(r) = V(°\r) + (— + —) vW(r) + 0{\) , (1) 

where (r) is the static potential, usually obtained by 
evaluating the expectation value of the Wilson loop. The 
static potential is well parameterized by the Coulomb 



plus linear term, 

y(°)( r ) 



+ or + (J, , 



(2) 



where a is the string tension and [i a constant [29j . On 
the other hand, the nonperturbatively expected form of 
V^(r) is not yet known, but leading-order perturba- 
tion theory yields V {1 \r) = ~C F C A a 2 J (4r 2 ) H, 13 
E3, where C F = 4/3 and C A = 3 are the Casimir 
charges of the fundamental and adjoint representations, 
respectively (beyond leading-order perturbation theory, 
see [3). 



II. PROCEDURES 

We work in Euclidean space in four dimensions on a 
hypercubic lattice with lattice volume V = L 3 T and lat- 
tice spacing a, where periodic boundary conditions are 
imposed in all directions. Writing the eigenstate of the 
pNRQCD Hamiltonian at O(m ) in the 3 ® 3* repre- 
sentation of color SU(3), which corresponds to the static 
quark-antiquark state, as |n) = |n; ri, 7^2} with the energy 
E n (r) [e.g., Eo(r) = (rjj^the spectral representation 
of (r) is expressed as 
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where g is the gauge coupling, E(fi) denotes the elec- 
tric field attached to the quark (i = 1) or the antiquark 
(i = 2), and AE n o = E n — Eo the energy gap. It is also 
possible to write Eq. © as the integral of the electric 
field strength correlator on the Wilson loop with respect 
to the relative temporal distance between two electric 
fields This is, in principle, measurable on the lat- 

tice, and the result is reduced to Eq. Q once the spec- 
tral decomposition is applied by using the transfer matrix 
theory, and the temporal size of the Wilson loop is taken 
to infinity [sfil ]. 

In our approach, the Polyakov loop correlation func- 
tion (PLCF, a pair of Polyakov loops P separated by a 



distance r) is adopted as the quark-antiquark source in- 
stead of the Wilson loop for the reason discussed below. 
Let us consider the field strength correlator on the PLCF, 

C(t) = ((ff 2 B(f i ,t 1 )-B(r' < ,ta)))c 

where the double brackets represent the ratio of expec- 
tation value ((• • •)) = (• • ■) P p»/(PP*), while (• • -)pp- im- 
plies that the electric field is connected to either of the 
Polyakov loop in a gauge invariant way. The relative tem- 
poral distance of two electric field operators is t = ti — 1\. 
The spectral decomposition of Eq. Q reads 
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where the last term represents terms involving expo- 
nential factors equal to or smaller than exp[— (AEio)T]. 
Thus, once Eq. (@J is evaluated via Monte Carlo simu- 
lations, we can determine the amplitude |(0|<7-E(ri)|n)| 2 
and the energy gap AE n g in Eq. (JSj by a fit and in- 
sert them into Eq. ©. It is easy to see that in the 
limit T — > oo we can write Eq. i n the integral form 
VW(r) = -(l/2)lim J -_ >oo / T dt tC(t), where r = t]T 
with arbitrary r\ <E (0,T/2]. 

The reason for using the PLCF is to compute Eq. Q 
with less systematic errors. The hyperbolic cosine in 
Eq. (JSJ is typical for the PLCF and we can control the ef- 
fect of the finite temporal lattice size on the field strength 
correlator automatically in the fit. Moreover, the error 
term of O(e~( A£;i0 ) T ) is already expected to be small for 
a reasonable size of T. By contrast, if one uses the Wil- 
son loop at this point, the spectral representation is just 
a multi-exponential function, and the leading error term 
is of O(e~' ABl0 ^ A *)), where At is the relative temporal 
distance between the spatial part of the Wilson loop and 
the field strength operator. Here, one cannot choose At 
as large as T, since the temporal extent of the Wilson 
loop is limited to T/2 because of the periodicity of the 
lattice volume. 

The only technical problem that arises when using the 
PLCF is how to obtain a signal for the field strength 
correlator in Eq. |@J , since the expectation value of the 
PLCF at zero temperature becomes exponentially small 
with increasing r, and the signal is easily washed out by 
statistical noise. In fact, it is almost impossible to ob- 
tain the signal of the PLCF at intermediate distances, 
say r ~ 0.5 fm, with the commonly used simulation al- 
gorithms. However, we find that this problem can be 
solved by applying the multi-level algorithm |20j with a 
certain modification as applied to the determination of 
the spin-dependent potentials 0,0] (see also for a 
similar application) . 

The basic procedure of the multi-level algorithm (re- 
stricted to the lowest level) is as follows. We first di- 
vide the lattice volume into several sublattices along the 



FIG. 1: Construction of the electric field strength correlator 
on the PLCF. Arrows at ri and r*2 represent the Polyakov 
lines for the static quark and antiquark. [■ ■ ■ ] denotes the 
sublattice average. 

time direction, where a sublattice consists of a certain 
number of time slices N ts \- The number of sublattices is 
N su b = T/iVtsi, which is assumed to be integer. In each 
sublattice we take averages of the components of the cor- 
relation function [components of the PLCF and of the 
field strength correlators, which are in the 3 ® 3* rep- 
resentation of SU(3)], by updating the gauge field with 
a mixture of heatbath (HB) and over-relaxation (OR) 
steps, while the spatial links on the boundary between 
sublattices remain intact during the update. We refer 
to this procedure as the internal update and denote the 
number of internal update as Mupd- Repeating the inter- 
nal update until we obtain stable signals for these com- 
ponents, we finally multiply these averaged components 
to complete the correlation function. Thereby the cor- 
relation function is obtained for one configuration. For 
a schematic understanding, see Fig. ^ which illustrates 
the computation of the electric field strength correlator 
on the PLCF. We then update the whole set of links 
without specifying any layers to obtain another indepen- 
dent gauge configuration and repeat the above sublat- 
tice averaging. Once N ts i and -/Vj up d are optimized for a 
given gauge coupling (3 and a maximal quark-antiquark 
distance of interest, the statistical fluctuations of observ- 
ables turn out to be quite small. Further technical details 
can be found in [Tgj ■ 



III. RESULTS 

Our simulations were carried out using the standard 
Wilson gauge action in SU(3) lattice gauge theory at 
j3 = 6.0 on the 20 4 lattice (the lattice spacing, determined 
from the Sommer scale ro = 0.5 fm, is a w 0.093 fm |20j). 
One Monte Carlo update consisted of 1 HB, followed by 
5 OR steps. For practical reasons (mainly to save com- 
puter memory) we set r = (r, 0,0). We employed the 
lattice field strength operator defined by ga 2 F flu (s) = 
[U^(s) — W {s)]/{2i) at the site s, where U^ v {s) are 
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FIG. 2: The electric field strength correlators on the PLCF 
at /3 = 6.0 on the 20 4 lattice for r/a = 5. The dotted lines 
are the fit curves with n max = 3 in Eq. 

plaquettc variables and constructed the electric field by 
ga 2 Ei(s) = ga 2 [F 4:i (s) + F<ti(s — i)]/2. In order to re- 
move self-energy contributions of the electric field we 
multiplied by the conventional Huntley-Michael factor, 
Ze^) which, however, only removes self-energy 

contributions at 0(g 2 ). This factor, which depends on 
r and also on the relative orientation of the electric field 
operator to f, was computed using the PLCF [l^. We 
obtained the value ZnAr) « 1.62. For a more precise 
value of Ze, see Ref. For the chosen value of -ZV ts i = 4 
we performed -/Vi up d = 7000 internal updates. Our total 
statistics was N con { = 60. 

In Fig. |21 we show the C(t) for the longitudinal and 
the transverse components, ((g 2 E x {fi,ti)E x (fi,t2))) c and 
({g 2 Ey(r i ,t 1 )Ey(f f i ,t 2 )}}c = ^E^r^t^E^nM)))^ re- 
spectively, where r/a — 5 is selected as an example. 
Note that the correlators are negative. Here, the sec- 
ond term of Eq. can be non-zero as the electric 
field is even under CP transformations. We computed 
{(gEi)) independently and found {(gE y )) = ((gE z )) = 0, 
while ((gE x )) ^ 0, which was then subtracted to ob- 
tain C(t). As it is impossible to determine the ampli- 
tudes and the energy gaps for all n > 1 with the lim- 
ited data points, we truncated the expansion in Eq. © 
at a certain n = n max . The validity of the trunca- 
tion was monitored by looking at x an d the stability 
of the resulting potential as a function of n max , where 
X was always defined with the full covariance matrix. 
We found that n max = 3 was optimal with the fit range 
t/a e [1,8] (equivalent to t/a G [12,19]). The system- 
atic effect caused by the truncation can be checked by 
simulating volumes with larger values of T and by in- 
creasing n max in the fit. However, from the experience 
of evaluating similar field strength correlators for the 
spin-dependent potentials [l9j ]. we expect that such an 
effect is already negligible compared to statistical errors, 
once three terms are included for T = 20 at (3 = 6.0. 



0.42 - 




2 4 6 8 



r/a 

FIG. 3: The potential at 0(l/m), V (1) (r). Dashed and solid 
lines are the fit curves corresponding to Eq. @ and Eq. J7|, 
respectively. 



Here, we employed two ways of the fit procedure; we 
fitted ({g 2 E x E x )) c and ({g 2 E y E y )) c separately, and fitted 
{{g 2 E-E)) c = {(g 2 E x E x )) c + 2((g 2 E v E y )) c simultaneously. 
The latter is based on the expectation that the energy 
gaps are the same for both correlators. We obtained 
xljN di = 1.1 for ({g 2 E x E x )) c and 3.0 for {{g 2 E y E y )) c , 
respectively, and the corresponding fit curves are plotted 
in Fig. Ndi is the number of degrees of freedom. The 
simultaneous fit yielded Xmm/^df = 2.2. In any case, 
the resulting potential was found to be the same within 
errors, which were estimated from the distribution of the 
jackknife sample of the fit parameters. For other dis- 
tances Xmin/-^df was smaller than in this example, and 
the results of the two fit procedures were consistent. 

We present the potential (r) in Fig. |3| where the 
result of the simultaneous fit is plotted. We see an in- 
creasing behavior as a function of r. We first tested 
whether this increasing behavior matches the expecta- 
tion from perturbation theory. Neglecting logarithmic 
corrections we fitted the data at r/a e [2, 5] to 

V^) = - C - + ^ (6) 

and found c' = 0.099(5) and aV = 0.401(1) with 
Xmin/-^Vdf = 6.6, where the fit curve is plotted in Fig. |3| 
(dashed line). Note that if we include the data at 
r/a = 6, x 2 becomes twice as large, while the fit pa- 
rameters are little affected. In order to check if this is a 
remnant of the perturbative behavior, we need data at 
smaller distances and perform a scaling test. At the mo- 
ment, what we can say is that the data at r/a > 5 are 
inconsistent with a pure 1/r 2 behavior. 

In trying to establish empirically the functional form 
of the r dependence, we employed several alternative fit 
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functions, and among them, we found that 

^M = -7+M w , (?) 

can describe the behavior of VW(r) reasonably well, 
where the coefficient c" has a dimension of mass. We 
took into account the data at r/a G [2,6] and obtained 
ac" = 0.081(4) and a 2 fx" = 0.417(1) with xLiJ N df = 
2.3, where the fit curve is plotted in Fig. [3] (solid line). 

As the potential V (1) (0 requires no matching coef- 
ficient |24l . l25j . in contrast to the spin-dependent po- 
tentials at 0(l/m 2 ), we can directly insert V^'(r) into 
Eq. Q and compare its relative magnitude with the 
static potential (r) for given quark and antiquark 
masses. For this purpose we may use the fit result of 
Eq. 0. By dividing V^ 2 ( r ) by the quark mass, where 
we set mi = rri2 = m for simplicity, we have a 1/r term 
with a dimcnsionlcss coefficient 2c" /m. For charmonium, 
m c = 1.3 GeV, we then find 2c"/m c = 0.26(1), which is 
93(5) % of the Coulombic coefficient of the static poten- 
tial, c = 0.281(5), in Eq. © 0. For bottomonium, 
m b = 4.7 GeV, we find 2c" /m b = 0.073(4), which is still 
26(2) % of c. It is certainly interesting to investigate the 
effect on heavy quarkonium spectroscopy. 

IV. SUMMARY 

We have investigated the relativistic correction to the 
static potential at 0(l/m) nonperturbatively by using 



lattice QCD Monte Carlo simulations for the first time. 
The key strategy here is to employ the multi-level algo- 
rithm for measuring the field strength correlator on the 
PLCF and to extract the potential by exploiting the spec- 
tral representation of the field strength correlators. This 
method allows us to obtain the potential with less statis- 
tical and systematic errors. The correction is found to be 
comparable to the Coulombic term of the static potential 
when applied to charmonium and to be one-fourth of the 
Coulombic term for bottomonium. 

Finally, we note that the field strength correlator 
obtained here can be used to compute one of the 
velocity-dependent potentials at 0(l/m 2 ), Vd(r), in the 
parametrization of Refs. since the spectral rep- 

resentation of Vd(r) consists of the same amplitudes and 
the energy gaps. We plan to present this result as well as 
the other velocity-dependent potentials at 0(l/m 2 ) in a 
separate publication. The first lattice result can be found 
in Ref. |2^ . 
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